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Abstract 

Spectral decomposition of dynamical equations using curl-eigenfunctions has been extensively used in 
fluid and plasma dynamics problems using their orthogonality and completeness properties for both linear 
and non-linear cases. Coefficients of such expansions are integrals over products of Bessel functions in 
problems involving cylindrical geometry. In this paper, certain identities involving products of two and 
three general solutions of Bessel's equation have been derived. Some of these identities have been useful 
in the study of Turner relaxation of annular magnetized plasma [S.K.H. Auluck, Phys. Plasmas, 16, 
122504, 2009], where quadratic integral quantities such as helicity and total energy were expressed as 
algebraic functions of the arbitrary constants of the general solution of Bessel's equation, allowing their 
determination by a minimization procedure. Identities involving products of three solutions enable 
expanding a product of two solutions in a Fourier-Bessel series of single Bessel functions facilitating 
transformation of partial differential equations representing non-linear dynamics problems into time 
evolution equations by eliminating spatial dependences. These identities are required in an ongoing 
investigation of fluctuation-driven coherent effects in a nonlinear djTiamical system. 
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I. Introduction 

Eigenfunctions of the curl operator, also known as Chandrasekhar-Kendall (C-K) functions [1] 
and their generahzation to other geometries, can be used as basis functions to expand vector fields in 
spectral decomposition methods because of their orthogonality and completeness as well as because of 
simplifications introduced in handling electromagnetic fields and fluid vorticity [2]. This is usually 
confined to linear equations because of the limitations involved in determining the spectral coefficients of 
such expansion. In cylindrical geometry, C-K functions have Bessel functions describing radial modes, 
sine and cosine functions describing axial modes in a periodic boundary condition and the azimuthal 
modes. If the djTiamical variable are expressed as Fourier-Bessel series in C-K functions (each C-K 
function representing an eigenmode), then their products in non-linear equations would be described as a 
double series in products of C-K functions. If a product of C-K functions itself can be expressed as 
Fourier Bessel series in C-K functions, then the spectral decomposition method can in principle be 
extended to non-linear equations as well, facilitating transformation of partial differential equations 
representing non-linear dynamics problems into time evolution equations by eliminating spatial 
dependences. This involves numerically evaluating and tabulating integrals involving products of three 
C-K functions. However, when large mode numbers, representing finer spatial scales, are involved, the 
integrands vary rapidly causing difficulties in numerical evaluation. Analytical expressions for such 
integrals would be very useful in dealing with large mode number situations enabling a deeper 
understanding of scaling relationships in the parameter space. If analytical methods are not available, a 
large computational effort required to generate a database of such integrals using arbitrary precision 
methods once for all would be quite useful. 

In this paper, an attempt is made to derive some integral identities involving products of Bessel 
fijnctions, encountered in simple problems with azimuthally symmetric geometry. Some of these 
identities have been useful in the study of Turner relaxation of annular magnetized plasma [3], where 
quadratic integral quantities such as helicity and total energy were expressed as algebraic functions of the 
arbitrary constants of the general solution of Bessel's equation, allowing their determination by a 
minimization procedure. A part of the following is included in the appendix of [3], which is re-derived 
below for the sake of completeness and continuity. 



II. Identities involving products of two functions: 

The fimctions 
Z„(x)^aJ„(x) + bY„(x) 



1 



where a and b are arbitrary constants, obey the following relations common to J„ and Y^^ , the Bessel 
fiinctions of first and second kinds with integral order n. 

-Zn(x) = Z„,(x) + Z„,,(x) 2 

X 

2— ^ = Z„_, X -Z„,, X 3 
ox 

Z„(x) = (-l)"Z„(x) 4 

If Up(ax) and Vq (Px) are two solutions of Bessel' s equation of order p and q, then the following identity 
is known [4] : 



(a^-P^)x-^^ 



Up(ax)V^(Px) dx 



= PxUp (ax) V^_, (Px) - axUp_, (ax) (Px) + (p - q)U^ (ax) (Px) 

Putting \J = Z^ (ctx); V = ^ (Px) , (where Z and Zcould in general have different values of a and b, 
a situation which might be encountered in some problems) and using 2 and 3 

Jxdx[z.(ax)^(Px)] = ^{pxZ.(„x)^..(Px)-c«Z..,(„x)^(Px)} 
- ^ Iz.(„x)W.^(p,).«.(<xx)l 



(a=-|3')l ' dx ^^'^ ' dx J 

In the case of a problem with annular cylindrical geometry [2], where the normalized inner radius 
corresponds to x=l, the normalized outer radius corresponds to x=A, (=outer radius /inner radius) and 
physical boundary conditions require that some of the fields (expressed as linear combination of Bessel 

fimctions for a given order n) vanish at the boundaries: Z^ (otA) = Z^ (a) = 2^ (PA) = 2^ (p) = Othe 

definite integral over the domain becomes 

f xdx[z. (ax)^ (Px)] = ' |(A=)M!!5>^.^!^!3?>^^ 

J I n\ j^yy^ jj i2a)[^ ^ dx dx dx dx J 

This orthogonality relation, which is a direct consequence of the boundary conditions, allows expressing 
quadratic integral forms, defined over an annular domain, involving physical fields, such as total energy, 
power dissipation, helicity, stress tensor in terms of the arbitrary coefficients a,b which define the linear 



combination in equation 1 and which can be determined at the end of the problem through an extremum 
principle. 

For subsequent discussion, Z and Z are taken to be identical. Substitute P=a+8 and expand in first power 
of 8 and take the limit as 8^0. 



jxdx[Z„(ax)]' =y [Z„(ax)]'-Z„^j(ax)Z„_j(ax) 
Forn=Oandn=l 
|xdx[ZQ(ax)J' 



l2 X 

T 



[Z„(ax)] +[z,(ax)] 



jxdxZj(ax)^ = [Zj(ax)]^ +[ZQ(ax)]^ - — Zg(ax)Zj(ax) 
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Multiplying 2 and 3 for x— >ax, recurrence relation is obtained for higher order integrals similar to 9, 10: 

12 2n 



Jxdx[z„^i(ax)] =|xdx[z„_i(ax)] -[Z„(ax)] 
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a 



Using relations — ZQ(ax) = -aZj(ax); — xZj(ax) = axZQ(ax) obtained from 2 and 3, the 
following can be derived 

|x%'(ax)dx 

= --x% (ax)Zo (ax) + jxP"M[xZi (ax)]' + jx" [Z^ (ax)]'dx 



12 



a 

For p=3 



|x^|[ZQ(ax)]' -[Zj(ax)]'|dx = a VZj(ax)Zo(ax)-a ^x^[Zj(ax)]' 
Using 8 with n=0 and 4, 



[Zo MT X = A|^[[Z, (ax)]V[Z, (ax) J 



13 



14 



Therefore 



x'[Zo(ax)] =x 



2 2 d X 



= X — 

dx 2 



Thus 



dx 2 

^2 



[Zo(ax)J+[Z,(ax)J 
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[ZQ(ax)] +[Zj(ax)] -x^ [ZQ(ax)] +[Zj(ax)] 



dx 



[Zo(ax)]'+[Zi(ax)]' 



16 



|x^[2[Z„(ax)]' + [Z,(ax)J 
Adding 16 and 13 

jx^[Zo(ax)]'dx = ^ [Zo(cQc)]'+[z,(ax)J +— x^Z^(ax)Zo(ax)-— x^[z^(ax )J 17 



3a' 



Subtracting 13 from 17 



jx^[Zj(ax)]^dx = ^ [Z(,(ax)J +[Zi(ax)]' - ^xX(ax)Zo(ax) + ^x^[Zj(ax)]^ 18 



3a' 



Multiplying 2 and 3 for x^ax, with a little manipulation, recurrence relation is obtained for higher order 
integrals similar to 17, 18 



[Z„,i (ax)]' dx = |x' [Z„_i (ax)]' dx 



4n x' 



+— ^ — 
a 



.2 , 2 



[Z„ (ax)] - Z„^i (ax)Z„_i (ax)J " [Z„ (ax)] 



2n 



19 



From 6, 



Z,(ax)Z„((3x) 



20 



= pxZj (ax)Z_j (px) - axZo (ax)Zo (px) + {ax)Z^ (px) = y;"^ (x) 
Therefore 

j xdYi"/ (x) = (a' - p' )| x'dxZi (ax) Z^ (px) - JdxZ; (ax) Z^ (px) 

= xYi'jP (x) + jpxZi (ax)Zi (px)dx + jaxZo (ax)Zo (px)dx- jz, (ax)Zo (px)dx 



21 



This yields the formula 



W,f (x) = jx^Zj (ax) Zq Ox)dx {a ^ p} 



2a' 

+ - 



22 



, , xZ^ (ax) Zp (px) - ^"^P xZp (ax) Z, (px) 

(a^-P^) (a^-P') 

The following result can be derived either from 22 or from 9 

1 2 

Jx^Zj(ax)Z(,(ax)dx = — x^[Zj(ax)J 23 

The identity 22 can be used to derive the following identities 

jx^dx[Z„(ax)Z„(px)] 

_ ax^Zi (ax) Zo (px) - px% (ax) (px) 

.^fMxZ„ (ax)Z. (px)-^^M^^ (a.)z„ (Px) 

|xMx[Z, (ax)Z, (Px)] = PxX(ax)Z,(px)-ax^Z,(ax)Z,(px) 

(a p j 

"^"^^ x'Z, (ax)Z, (px) - ^^"^^ x% (ax)Zo (px) 25 



4p(a'+p') ^ . ^ 4a(a'+p') x 

- , ; ,y (Px)Z„ (ax)+ \ ^ xZ, (ax)Z„ (px) 

Some of these results have been used in ref [3]. 

III. Identities involving products of three functions: 

In the following, functions representing indefinite integrals will be written with explicit ilinctional 
dependence (x). Define 



C (x) - |xZ„ (ax)Zo (Px)Z„ (Yx)dx 



26 



irf/ (x) - JxZ, (ax)Z, (Px) (Yx)dx 
C (x) - JxZ, (ax)Z, (px)Z, (Yx)dx 
C (x) - jxZo (ax)Zo (Px)Z, (Yx)dx 
Kit (x) - jz, (ax)Z, (px)Zo (Yx)dx 
Kt (x) ^ jz, (ax)Z, (Px)Z, (yx)dx 

W«P (x) = f xdx[Zi (ax)Zj (px)] = / .J PxZ, (ax)Zo (px) - oxZ^ (ax)Z, (px)} 

(a -p j 

W;oP (x) - jxdxZo (ax)Zo (Px) = ^-^l-^{axZ, (ax)Zo (px) - PxZ^ (ax)Z, (px)} 
From 26 

C (x) - jz„ (ax)Z„ (px)xZ„ (yx)dx = ^jz^ (ax)Z„ (px)d[xZ, (yx)] 
= ^[xZi (yx)Zo (ax)Zo (px)] - ^ jxZ^ (Yx)d[Zo (ax)Zo (px)] 
= \xL, (yx)Z, (ax)Z„ (Px)] A\t (x)+ ^ if- (x) 



27 
28 
29 
30 
31 

32 
33 



34 



Therefore, 

allt (x) + pif- (x) = yl2 (x) - [xZ, (yx)Z„ (ax)Z„ (Px)] 

By cyclic substitution a^P, P^y, y^a two additional equations are obtained, which can be 
written in the matrix form 








Y 


"C(x)" 




a 


Y 





C(x) 







P 


a 


C^x) 





aC(x)-[xZ,(ax)Z„(Px)Z„(yx)]" 
PC (x) - [xZ, (Px) Z„ (yx) Z„ (ax)] 
yC (x) - [xZ, (yx)Zo (ax)Zo (px)] 



35 



36 



This has the solution 



C(x) = Iot(x)(2apr(a^+P^-y^) 

+ (2ap) ' [yxZj {y^)Z^ (ax)Zo (px) - axZ^ (ax)ZQ (px)Zo (yx) - pxZj (px)Zo (yx)Zo (ax)] 
Two similar expressions can be written by cyclic substitution. 
For the special condition (a^ + — ) = , 37 becomes 

C(x){y^V^} 

= (2ap)"' [yxZi (yx)Zo (ax)Zo (px)-axZi (ax)Zo (px)Zo (yx)-pxZi (px)Zo (yx)Zo (ax)] 
An alternate (but not independent) relation between the two integrals can be obtained as follows 
C(x) = |xZ„(ax)Z„(px)Z„(yx)dx = |z„(yx)d[wr(x)] 

= (x)Z, (yx) - j W^^ (x)d[Z„ (yx)] 39 

= (x)Z„ (yx) + ^^C^ (x) - lf[o (x) 

(a -p j I" ~P j 

Iqqq (x) satisfies the following relationship 

A differential relationship can be demonstrated between the two integrals I^J (x)and V^f^ (x) : 
^C(x) = -i£aC(x);Ac(,)._iApC„ 

This can be used to eliminate one of the integrals giving a differential equation for (x) 



^C(x)-F(a,p,y)C(x) = G(a,p,y,x) 

F(a.p.y)= ^°(f 7'-°') =A{_lLog[(a-.p--rr-4aV 
G(a,|3,Y,x)= ^ 



y(a^+P^-y^)x%(ax)Z„(Px)Z,(yx) 
+ (a^-P^-y^)x%(ax)z4Px)Z„(yx) 
-(a^+p^-y^)x^Zj(ax)Z^(px)Z„(yx) 
-2aPyx^ZQ (ax)Zj (px)Zj (yx) 
+2a'px'Zi (ax)Zi (Px)Zo (yx) 

Although the solution of 42 can be formally written down, no simple formula is found for (x) . 

Starting from 28, 

C(^) = ^jz,(ax)Z,(Px)d[xZ,(yx)] 

1 OC r 
= -xZj(yx)Zj(ax)Zj(Px) xZj(yx)Zj(Px)dxZQ(ax) 

y y 

-— [xZj(yx)Zj(ax)Zo(Px)dx + — [Zj (yx)Zi (Px)Zi (ax)dx 

y J y J 

Therefore 

yC (x) + aC" (x) + PIJ^P (x) = xZ, (yx)Z, (ax)Z, (Px) + 2T^^ (x) 
Using 37 

^aP. (^) = (^'-Pj-y') xZ, (ax)Z„ (px)Z„ (yx) + ^-^^^xZ, (Px)Z, (yx)Z, (ax) 
4py 4ya 

+ ^ ^ xZ, (yx)Z„ (ax)Zo (Px) - - xZ, (yx)Z, (ax)Z, (px) 



-C(x)- 



(a'+p'-y')'-4ay 
4aPy 



( 2 2 2 2 2 

a +P -y ) — 4a P =0,45 becomes 



Kit/(x){y^|a±p|} 

= ^""^ ^ xZ, (ax) Z„ (Px)Z, (yx) + '^^^ ^ xZ, (Px)Z, (yx)Z„ (ax) 46 

^ (y -a -(3 ) ^^(y^)z„(ax)Z„(px)-lxZ,(Yx)Z,(ax)Z,(px) 
4ap 2 

In a similar maimer, the following identities can be proved 

(x) + aC, (x) + pi^^f (x) = -xZ, (ax)Z, (px)Z„ (yx) + K,t 47 



I,, (x)-W„ (x)Z.(Yx)-p3^V +pZ^Ioo> ^ 



Pya 
10 



48 



PC -KrfJ - yC (x) = Zo (yx)xZ, (ax)Z. (Px) 49 



Cyclic substitution in 49 leads to three simultaneous equations which can be used to express Iq^ J etc in 
terms of K"fj and I"f/ and their cychc variations. 



2api^P/ (x) = aKf- (x) + pK-^ (x) - yK^f; (x) + I^f/ (x)(a^ + p^ - y^ ) 

+axZQ (ax)Zj (Px)Zj (yx) + PxZ^ (Px)Zj (yx)Zj (ax) - yxZ^ (yx)Zj (ax)Zj (Px) 



IV. Case of Bessel functions of first kind 

When the Z functions in the previous section are chosen as Bessel fimctions of the first kind, 
definite integrals corresponding to equations 26-33 within the limits {0,1} can be defined. Additionally, 
a, P, y are assumed to be zeros of Bessel functions: cx, ^ p which is the p* zero of Bessel function of 

first kind of order q (q=0 or 1). Define 

,Cr-ixJ.(vx)j.(vx)j,(j^,px)dx 51 





From 37 



mnp mnp 

q'^llO ~ q'^OOO 



From 45 



^ Jq,mJq,n 



53 



pvmnp pi mnp ^ 



^Jq,mJq,nJq,p 



54 



No formula has yet been found for ■ An approximate scaling relations for j C is obtained in the 



next section. 



The importance of these integrals may be appreciated from the following. The product of two 
Bessel fiinctions may be expanded into a Fourier-Bessel series: 



jj(i„x)Jk(ipx)=|;c.7J,(j,px) 



55 



The coefficients of this expansion are given by 

c,7=2[ji4lp)]>7 56 
As an illustration, consider the generalized Ohm's law for a two-fluid plasma with electron inertia [5] 



SoCOp I di 



K + VXD = I] J -^ -jjxr)— vpg — V •■; r- f> 

Using Maxwell's equations, this can be written as 



2T - ----- 

— ^ + v — + co„K + co>xr> + vc VxB + c^Vx Vxli- 



-I- V -^JV -I- VJ 



'1 



57 



58 



e - - _ - - 
J xB + 8„'V-|Jv-i- vj 



Some 



ens„ 



V-(^JJj + 



J(^J -Vnj 



Expressing the electric and magnetic fields in terms of C-K functions in cylindrical geometry with 
azimuthal sjmimetry having time dependent mode amplitudes, a system of time evolution equations can 
be obtained whose left hand side consists of a linear differential operator acting on the mode amplitudes 
and right hand side consists of a "sum over modes" of terms involving products of quantities related to 
the mode amplitudes weighted by coefficients of the form given by 51 and 52. In the case being 

considered, it turns out that all the coefficients are of one of the three form: ^^'^w 1^1™'' or iDj"J7 . 



This procedure allows taking into account effects due to the electron momentum convection term 
V ^^JJj which is of order c^ai^i as compared with the Hall effect term, (, being a gradient scale 

length and cOp is the electron plasma frequency. This term is usually too stiff to be used in the 

conventional approach using partial differential equations. The difficulty in niunerical treatment of this 
term in a system of difference equation is converted into the difficulty of evaluating the coefficients of the 

form iC™o'' or jC,™''or jD^^for large mode numbers (finer spatial scale). Identities 53 and 54 reduce 
the problem to evaluating only jC^qq'' as a function of mode numbers. The elecfron momentum convection 
term, being a gradient, should give a higher contribution at higher mode numbers. Since the sum over 
modes has to be evaluated for a very large number of modes (-10^), a relationship showing how ,CJJU5f 

scales with the mode numbers m, n, p is very desirable to enable estimation of the net effect of the 
electron momentum convection term and similar other terms to observable average effects. 

V. Approximate scaling relation for ,c^f 



mnp 
000 



In the absence of an analytic formula, an approximate scaling relation is derived below for jC 

for the case of large mode numbers, where the integrand oscillates very rapidly creating difficulty with 
numerical evaluation of integrals. For smaller mode numbers, the coefficients can be numerically 
evaluated and tabulated once for all. 

The starting point is the result that p* zero of Bessel function of first kind of order 1 is jp~p7i for large p 
and the asymptotic approximation for Bessel fiinctions of first kind for integral orders [6] for large values 
of the argument: 



J (x)D cosfx-^nTi-jTi:)] for large X 59 

V Tlx ^ 



This leads to the following useful approximate formula 



3pp-j^%(jp'^)jo(jpx)d^= 





2 



(-ir" (-1) 



p+p 



(jp'-jp) (jp'-jp) (Jp'+jp), 
1 



for p^p' 



for p = p' 



The triple product integral can be approximated as follows: 



C0s[7i(m + n + p)^ - ^(i + j + k)7i - Jti] 
^1 j. +cos[7i(m-n-p)^-^(i-j-k)7i + |7i] 



4 Ti-'^nmp i ^ 



It can be shown that 



|-^cos(Prt^ + Q7r) 
-v^ 



+ C0s[7i(m + n-p)^-^(i + j-k)7i:-^7i] 
+ cos[7i(m-n + p)^- j(i- j + k)7i-^7i] 



= 2cos(Q7r) 



5p,o + (l-5p,o)- 



c(^ 



-2sin(Q7r)sign(P) 



where the Fresnel Integrals are defined by 

t t 
S(t) = |sin(7it" / 2)dt'; C(t) = jcos(7it" / 2)dt' 



Applying this to 6 1 



^ijk ~ T 3 r 



mnp 



(_i(i + j + k)7l-|7l)^ 

^|(m + n + p)| 
sin(-i(i + j + k)7i;-|7r)s|^(m + n + p) j 



c(^2|(m + n + p) 



J(m + n + p)| 



+ COs(-Y(i- j-k)7i + ^7T:) 



^(m-n-p),0 + (-^ ^(m-n-p),0 ) 



C(^2|(m-n-p)|) 



-sign(m-n-p) 



J(m-n-p)| 
sin(-i(i-j-k)7i + i7r)s|^2|(m-n-p)|j 



^|(m-n-p)| 



+ COs(-^(i + j-k)7t-^7l) 



^(m+n-p),0 



CL 2 



(m + n - p) j 



-sign(m + n - p) 



J(m + n-p)| 
sin(-i(i + j-k)7i-i7i)s|^2|(m + n-p)|j 



-cos(-^(i-j + k)7r->) 



J(m + n-p)| 

^(m-n+p),0 ~ ^{m-n+p),0 ) " 



cU2|(m-n + p)| 



63 



-sigii(m-n + p) 



J(m-n + p)| 

sin(-i(i- j + k)7i-|7r)s|^2|(m-n + p)|j 
J(m-n + p)| 



VI. Numerical Validation 

The identities derived above have been validated by computing the left and right hand sides separately for 
varying arguments. The difference in the LHS and RHS ranges from 10"** to 10"^^; the poorer 
correspondence is correlated with the divergence of the Bessel function of second kind at small values of 
the argument. The approximate nature of relation 60 and numerical evaluation of the derivative in 41 limit 
the agreement to v^^ithin 10"^-10 The most difficult challenge is posed by validation of 63, for the case of 

jC^'' . The approximation represented by 59 is expected to be more accurate for larger mode nimibers 
while the numerical evaluation of C™'' becomes very difficult for large mode numbers. As a 
compromise, numerically integrated values and those calculated from the approximation formula are 



compared in Table 1 and fig. 1 for various combinations of the mode indices. Above mode number 150, 
the numerical integration converges very slowly and for some mode combinations, fails the convergence 
test. The agreement between the numerical integration and formula 63 is seen to be better at higher mode 
numbers than at lower mode numbers, as expected. However, even at higher mode numbers, there is some 
disagreement, when the numerical integration gives a negative number while the formula doesn't. It is not 
known whether this is due to numerical error in the integration because of a highly oscillatory integrand 
or because of the error in the approximation. The inclusion of higher order terms in 59 is not expected to 
improve matters very much, firstly because their analytical integration in 61 is not straightforward and 
secondly because their contribution at higher mode numbers is expected to diminish faster than that of the 
dominant term. 
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Table 1: 



TABLE-1 


m 


n 


P 


LHS 


RHS 


44 


23 


63 


4.557E-05 


3.140E-05 


20 


20 


20 


9.061 E-05 


8.071 E-05 


22 


22 


22 


7.508E-05 


6.689E-05 


25 


25 


30 


5.265E-05 


4.659E-05 


30 


30 


30 


4.065E-05 


3.627E-05 


22 


89 


31 


-7.053E-10 


9.828E-06 


40 


45 


50 


1.866E-05 


1.787E-05 


37 


77 


57 


1 .590E-05 


1 .408E-05 


47 


61 


87 


1.162E-05 


1 .032E-05 


29 


47 


57 


2.342E-05 


2.158E-05 


40 


40 


40 


2.297E-05 


2.053E-05 


50 


50 


50 


1.474E-05 


1.320E-05 


60 


60 


60 


1.025E-05 


9.195E-06 


70 


70 


70 


7.543E-06 


6.770E-06 


80 


80 


80 


5.781 E-06 


5.195E-06 


90 


90 


90 


4.571 E-06 


4.112E-06 


100 


100 


100 


3.704E-06 


3.335E-06 


105 


85 


97 


4.139E-06 


3.899E-06 


100 


50 


20 


-1.114E-09 


7.515E-06 


120 


120 


120 


2.575E-06 


2.321 E-06 


130 


130 


130 


2.195E-06 


1 .980E-06 


140 


140 


140 


1.893E-06 


1.708E-06 


150 


150 


150 


1.649E-06 


1.489E-06 


160 


160 


160 


1.450E-06 


1.310E-06 


160 


80 


70 


-1.052E-08 


2.131E-06 


170 


170 


170 


1 .285E-06 


1.1 61 E-06 


200 


200 


185 


9.807E-07 


9.184E-07 


200 


200 


200 


9.286E-07 


8.401 E-07 


200 


100 


185 


1 .749E-06 


1 .582E-06 



Comparison of LHS and RHS of 63. The LHS is 
calculated by numerical integration over the product of 
three Bessel functions. RHS is the approximation 
formula. The trial values of mode numbers m,n,p are 
chosen so as to get a feel of how the integral behaves in 
various situations. 
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Fig-1 

Comparison of the approximation formula witli tlie numerical integration. The trend line has 
slope 1. The approximation formula tends to underestimate the integral. The three points nearly 
on the y-axis actually have a value close to zero and the numerical integration shows errors in 
convergence. 



